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Localized plasmons in metallic nanostructures have been widely used to enhance nonlinear op¬ 
tical effects due to their ability to concentrate and enhance light down to extreme-subwavelength 
scales. As alternatives to noble metal nanoparticles, graphene nanostructures can host long-lived 
plasmons that efficiently couple to light and are actively tunable via electrical doping. Here we 
show that doped graphene nanoislands present unique opportunities for enhancing nonlinear optical 
wave-mixing processes between two externally applied optical fields at the nanoscale. These small 
islands can support pronounced plasmons at multiple frequencies, resulting in extraordinarily high 
wave-mixing susceptibilities when one or more of the input or output frequencies coincide with a 
plasmon resonance. By varying the doping charge density in a nanoisland with a fixed geometry, 
enhanced wave mixing can be realized over a wide spectral range in the visible and near infrared. 

We concentrate in particular on second- and third-order processes, including sum and difference fre¬ 
quency generation, as well as on four-wave mixing. Our calculations for armchair graphene triangles 
composed of up to several hundred carbon atoms display large wave mixing polarizabilities com¬ 
pared with metal nanoparticles of similar lateral size, thus supporting nanographene as an excellent 


material for tunable nonlinear optical nanodevices. 
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I. INTRODUCTION 

Nonlinear optical phenomena arise from the coupling 
between two or more photons, mediated by their interac¬ 
tion with matter, to produce a photon with a frequency 
that is a linear combination of the original photon fre¬ 
quencies. These processes are responsible for many sig¬ 
nificant advances in laser-based technologies mm , most 
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of which rely on phase-matching of intense electromag¬ 
netic fields in extended bulk crystalline media to enable 
efficient frequency conversion. Now, as the accessibility 
of nanostructured materials offered by modern nanofab¬ 
rication techniques continues to increase, so does the in¬ 
terest in mastering nonlinear optics on subwavelength 
scales. Indeed, nonlinear optical nanomaterials find di¬ 
verse applications, including optical microscopy 
biological imaging/detection [5l[6], and signal conversion 
in nanoscale photonic devices [7Hin]. 

Inherently, nanostructured materials possess small vol¬ 
umes that limit their interaction with optical fields. For¬ 
tunately, this can be compensated by high oscillator 
strengths provided by either quantum confinement ef¬ 
fects m or the intense near-fields generated by localized 
surface plasmons [3 ng. In particular, the near-field 
enhancement associated with plasmons in noble metal 
nanostructures has been widely used with the purpose of 
enhancing the nonlinear response of surrounding dielec¬ 
tric materials Eiiiniiia. Plasmons also display strong 
intrinsic optical nonlinearities, clearly observed in metal 
nanoparticles nMs]. 

The plasmonic response of a metal nanostructure can 
be tailored by its size, shape, and surrounding environ¬ 
ment m, as well as by combining two or more interact¬ 
ing nanostructures. In this manner, a metallic nanocom¬ 
posite can be engineered to possess multiple resonance 
frequencies, enabling the optimization of nonlinear wave 
mixing between multifrequency optical fields dzi [nj. 
Unfortunately, plasmons are hardly tunable in metals af¬ 
ter a structure has been fabricated [19], thus severely 
limiting the choice of frequency combinations. 

As an alternative plasmonic material, electrically 
doped graphene has been found to support long-lived 
plasmonic excitations that efficiently couple to light and 
are actively tunable by changing the density of charge 
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FIG. 1: Plasmon-enhanced wave mixing of coincident pulses, (a) Illustration of an armchair triangular graphene 
nanoisland containing N = 1260 carbon atoms (side length of 8.4 nm) and doped with charge Q = 3 e, illuminated by two 
collinear light pulses of central energies hui = 0.41 eV and huj2 = 0.66 eV, each having 200 fs FWHM duration and peak 
intensity 10^^ W/m^. (b) Spectral density of the induced dipole moment under excitation by light polarized along a direction 
parallel (upper panel) or perpendicular (middle panel) to one of the nanotriangle sides (dipole along the incident electric-field 
direction). The solid black curves correspond to the spectra obtained upon excitation by coincident pulses, while the filled 
curves show the spectra produced by individual pulses of central frequency uji (red) or uj2 (green). The linear absorption 
cross-section of the nanoisland (normalized to its area) is presented in the lower panel (solid curve), where the arrows indicate 
the central frequencies of the exciting pulses and we compare it with the absorption of the undoped island (dashed curve). 


carriers [2DH32]. Strong intrinsic nonlinearities have 
also been observed in this material [33H35] . which could 
be further enhanced by plasmons [Ml EZl- Recently, 
plasmon-assisted second-harmonic generation (SHG) and 
down conversion with good efficiencies at the few-photon 
level have been shown to be possible when the funda¬ 
mental and the second-harmonic are simultaneously res¬ 
onant with plasmons in a graphene nanoisland [38]. We 
have predicted that a similar mechanism can lead to un¬ 
precedentedly intense SHG and third-harmonic genera¬ 
tion (THG) in nanographene [39]. Although experimen¬ 
tal studies have so far demonstrated strong plasmons 
at mid-infrared and terahertz frequencies PUHSU] . fur¬ 
ther reduction in the size of the islands down to less 
than 10 nm should allow us to reach the visible and 
near-infrared (vis-NIR) regimes [32l . In particu¬ 

lar, commercially-available polycyclic aromatic molecules 
sustain plasmon-like resonances that are switched on and 
off by changing their charge state m- Additionally, re¬ 
cent progress in the chemical synthesis of nanographene 
[33H37] provides further stimulus for the use of this ma¬ 
terial to produce electrically tunable vis-NIR plasmons, 
as well as their application to nonlinear optics at the 
nanoscale. 


In this work, we investigate nonlinear optical wave 
mixing in doped nanographene. Specifically, we pro¬ 
pose a scheme for optimizing wave mixing among mul¬ 
tifrequency optical fields that utilizes the plasmons sup¬ 
ported by a doped graphene nanoisland. Using a tight- 
binding description for the electronic structure combined 
with density-matrix quantum-mechanical simulations, we 
demonstrate that efficient wave mixing is achieved in an 
island when the incident and/or the mixed frequencies 
are coupled with one or more of its plasmons. By ac¬ 
tively tuning the doping level in a graphene nanoisland 
of fixed geometry, a wide range of plasmon-enhanced in¬ 
put/output mixing frequency combinations can be real¬ 
ized. 


II. RESULTS AND DISCUSSION 

We study wave mixing of coincident light pulses in 
Fig.[T] for a triangular armchair-edged graphene nanois¬ 
land containing N = 1260 carbon atoms (side length pf 
8.4 nm) and doped with three electrons (doping density 
9.1 X 10^^ cm“^, corresponding to an equivalent extended 
graphene Fermi energy Ep = 0.35 eV). Additionally, we 
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FIG. 2: Sum frequency generation We show the nonlinear polarizability + 002 )^ corresponding to sum-frequency 

generation (SFG), for an armchair-edged triangular graphene nanoisland containing N = 330 carbon atoms and doped with (a) 
three and (b) six additional electrons. The upper plots show the linear absorption cross-section, while the lower panels show 
the polarizabilities as the incident field frequencies uji and UJ 2 are varied. 


assume a conservative inelastic lifetime r = 33 fs (i.e., 
linewidth hT~^ = 20meV). From the linear absorption 
cross-section of the nanoisland, we identify several promi¬ 
nent plasmons for incident photon energies. In particu¬ 
lar, we concentrate on huji = 0.41 eV and fiuj 2 = 0.66 eV. 
The first of these plasmons is highly tunable upon elec¬ 
trical doping, as it is switched on when moving from 
Q = 0 to Q 7 ^ 0 charge states (see lower panel in Fig. 

and its frequency increases with Q (see Fig. in 
the Appendix). The plasmon at UJ 2 is comparatively less 
tunable. Plasmon-enhanced optical wave mixing in the 
nanoisland is demonstrated in Fig. Ef- where the spectral 
decomposition of the induced dipole moment is shown 
for excitation by collinear Gaussian pulses with central 
frequencies uji and 002 - We show the response for inci¬ 
dent polarization aligned with an edge of the nanoisland 
(upper panels) or perpendicular to an edge (middle pan¬ 
els), and in each case the dipole is calculated along the 
direction of polarization. Further, in Fig. we super¬ 
impose the spectra obtained from coincident pulses with 
the spectra for excitation of the nanoisland by each of 
the individual pulses in isolation, which exhibit polar¬ 
ization features that oscillate at harmonics of the funda¬ 
mental frequencies {nuji and nco’ 2 ). The dual-pulse spec¬ 
trum exhibits these features in addition to polarization 
components produced via sum and difference frequency 
generation {uji ± 002 ) and degenerate four-wave mixing 
{uji ± 2 ci ;2 or 2uji ± co’ 2 ), along with various other com¬ 
binations of harmonic generation and wave mixing. We 


note that for polarization along the edge of a nanotrian¬ 
gle, inversion symmetry in this direction prevents even- 
ordered nonlinear processes from occurring, such as SHG 
and sum/difference frequency generation. Incidentally, 
mixing of less tunable plasmons (e.g., huj 2 with the reso¬ 
nance at 1.20 eV) produces less intense nonlinear features 
(see Fig. in the Appendix). 

While the response of a graphene nanoisland to ul¬ 
trashort pulses can provide information on the relative 
strengths of the nonlinear processes for excitation at spe¬ 
cific frequencies, a quantitative analysis of each nonlin¬ 
ear process, along with its optimal plasmonic enhance¬ 
ment, is best provided by studying the response under 
continuous-wave (cw) illumination. In Fig. we con¬ 
sider sum-frequency generation (SFG) in a nanoisland 
containing N = 330 atoms (side length of 4.4 nm), for 
doping with either three (Fig. If) or six (Fig. [^) addi¬ 
tional charge carriers. The upper panels in Fig. [^show 
the linear response of the nanoisland for the two dop¬ 
ing levels considered, which show multiple plasmon res¬ 
onance peaks at low doping that converge to a single, 
stronger feature as the doping level increases. The SFG 
polarizability 002 ) is presented as a function of 

the two applied field frequencies, enabling exploration of 
all possible frequency combinations that may result in 
plasmon-enhanced wave mixing. 

The input frequencies at which SFG is enhanced are 
found to coincide with the plasmons, corresponding to 
the horizontally and vertically aligned features in the 
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FIG. 3: Four-wave mixing We show the nonlinear polarizability a^^\2(jJi — (JJ 2 ), corresponding to four-wave mixing, for the 
same nanoisland and doping conditions considered in Fig. 


contour plot of 002 )\ in Fig. Strong en¬ 

hancement is observed where these features intersect, 
driven by plasmons excited at both fundamental frequen¬ 
cies. Prominent features also follow frequencies satisfy¬ 
ing UJ 2 ojI = where uOp denotes one of the plasmon 
frequencies of the nanoisland, indicating plasmonic en¬ 
hancement at the output (sum) frequency. Finally, since 
SHG can be considered as a special case of SFG, we also 
note plasmonic enhancement of SHG for frequencies sat¬ 
isfying UJi = UJ 2 = UJp. 

Similar conclusions are drawn from frequency- 
difference generation (i.e., — co’ 2 )), where we ob¬ 

serve enhancement when either one of the incident fre¬ 
quencies or the output frequency resonates with a plas¬ 
mon (see Fig. in the Appendix). 

Plasmonic enhancement of four-wave mixing is inves¬ 
tigated in Fig. for the same graphene nanoisland and 
doping conditions considered for Fig. with an output 
frequency 2uji — Co’ 2 . We find a similar enhancement in 
a^^\2uji — Co’ 2 ) when either of the fundamental frequen¬ 
cies of the incident fields are resonant with plasmons in 
the nanoisland, although in this case the enhancement 
favors co’i = co’p, as co’i is involved twice in these specific 
wave mixing processes. Here we find increased four-wave 
mixing following the frequencies 002 = 2co’i — cjp, once 
again indicating enhancement at the output frequency. 
For Co’i = ^2 we have fully-degenerate four-wave mixing, 
where we are actually investigating the third-order polar¬ 
izability (cji) contributing to the linear response (i.e., 
the Kerr effect). In this case, frequencies satisfying the 
condition cji = ci ;2 = (^p produce extremely high polar¬ 


izabilities, as the three input frequencies and the output 
frequency are all equally amplified by the same plasmon 
resonance. We have also examined four-wave mixing with 
output frequency 2co’i +Co’ 2 , leading to similar conclusions 
on plasmon-assisted enhancement (see Fig. in the Ap¬ 
pendix). 

As noble metal nanoparticles possess the highest non¬ 
linear polarizabilities per atom measured experimentally 
m, we compare them with the above results for doped 
nanographene nonlinear polarizabilities. In particular, 
we contrast data available in the literature for metal 
third-order polarizabilities (uj) with those calculated 
here for graphene. In Fig. & we show the maximum 
polarizability |(a^^^(co’)| for graphene nanoislands of in¬ 
creasing size, and for a fixed doping ratio of one electron 
per every hundred carbon atoms, along with experimen¬ 
tal data obtained from degenerate four-wave mixing ex¬ 
periments on gold, silver, and copper nanoparticles of 
comparable or greater sizes. Even through metal parti¬ 
cles have much larger volumes than graphene nanoislands 
for a given lateral size L, we find the latter to exhibit 
much larger nonlinear response, reaching two orders of 
magnitude higher values. Obviously, we are comparing 
theory for graphene with experiments for metals, so this 
is a preliminary conclusion, which should be examined 
experimentally. Nonetheless, we are using a conserva¬ 
tive inelastic lifetime for nanographene in our calcula¬ 
tions (r = 33 fs), whereas much longer lifetimes could 
be actually encountered in practice, based upon recently 
measured plasmon resonances [26l [52] , and it should be 
noted that the nonlinear polarizability of order n scales 
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FIG. 4: Third-order response of nanographene compared to that of noble metal nanoparticles, (a) The maximum 
nonlinear polarizabilities of graphene nanoislands (solid triangles) as a function of their side length L are compared with 

measured values of noble metal nanoparticles of similar diameter reported in the literature. The color scale indicates the incident 
photon energy at which the maximum graphene nonlinear polarizabilities are found. The semi-transparent solid lines are guides 
to the eye, added to indicate the scaling with L. (b) Linear absorption spectra for the armchair-edged graphene nanotriangles 
considered in (a), distinguished by the number of hexagons along the edge rih , for which a fixed doping density of one electron 
per 100 carbon atoms is maintained. The frequencies at which the maximum values of (a^^^(cj) occur are indicated by the 
solid triangles, which follow the color scale in (a). The experimental data used in (a) is extracted from Refs. BBHSU, where 
the polarizabilities are obtained upon multiplication of the reported third-order susceptibility by the particle volume. 

Specifically, we show data obtained via fully degenerate four-wave mixing measurements of several noble metals performed 
at different wavelengths: gold nanoparticles at ~ 530nm [IH] (hollow squares), silver nanoparticles at ~ 420nm |49j (hollow 
diamonds) and near ~ 410nm [51] (solid diamonds), and copper nanoparticles at ~ 570nm [iO] (hollow circles) and ~ 530nm 
[50] (solid circles). 


roughly as [39] . The linear absorption spectra for the 
nanoislands that we compare with metal nanoparticles 
are presented in Fig. if, illustrating how the plasmon 
resonances vary as the nanoisland size increases while a 
fixed doping density is maintained. 


III. CONCLUSIONS 

In conclusion, graphene nanoislands can be used to 
realize nonlinear wave mixing on the nanoscale with ex¬ 
traordinarily high efficiencies, surpassing those of metal 
nanoparticles of similar lateral sizes. The large magni¬ 
tudes of the predicted nonlinear polarizabilities are at¬ 
tributed to plasmonic enhancement in these nanostruc¬ 
tures. Wave mixing is further enhanced by simultane¬ 
ously exploiting multiple plasmonic resonances, enabling 
two incident fields with distinct frequencies to indepen¬ 
dently couple to plasmons, and additionally tuning the 
mixed frequency to yet another plasmon. These plas¬ 
mons can be tuned by changing the number of dop¬ 
ing charge carriers, adding another advantage with re¬ 
spect to conventional plasmonic metal response. Inter¬ 


estingly, by considering islands formed by a few hun¬ 
dreds or thousands of carbon atoms, we predict tunable 
plasmonic response and plasmon-induced enhanced non- 
linearities within the visible and near-infrared spectral 
ranges. Our results configure a new platform for the de¬ 
velopment of nanoscale nonlinear optical devices based 
upon graphene nanostructures with lateral dimensions of 
only a few nanometers, such as those that are currently 
produced by chemical synthesis [i5M7] . 


IV. METHODS 

We describe the low-energy (< 3eV) optical response 
of graphene nanoislands within a density-matrix ap¬ 
proach, using a tight-binding model for the 7r-band 
electronic structure ISSl IM]- One-electron states \ipj) 
are obtained by assuming a single p orbital per car¬ 
bon site, oriented perpendicular to the graphene plane, 
with a hopping energy of 2.8 eV between nearest neigh¬ 
bors. In the spirit of the mean-field approximation m, 
a single-particle density matrix is constructed as p = 
^jj' Pjj'\^j){^j'\^ Pjj' time-dependent com- 
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plex numbers. An incoherent Fermi-Dirac distribution of 
occupation fractions fj is assumed in the unperturbed 
state, characterized by a density matrix pPj, = Sjj/fj^ 
whereas the time evolution under external illumination 
is governed by the equation of motion 

| = ( 1 ) 

The last term of Eq. [^describes inelastic losses at a phe¬ 
nomenological decay rate 1 /r. We set hT~^ = 20meV 
throughout this work (i.e., r = 33 fs), corresponding 
to a conservative Drude-model graphene mobility /i 
460cm^V“^s“^ for a characteristic doping carrier den¬ 
sity 4 X 10^^ cm“^ (i.e., one charge carrier per every 100 
carbon atoms). The system Hamiltonian H = Htb — 
consists of the tight-binding part i^xB (he., nearest- 
neighbors hopping) and the interaction with the self- 
consistent electric potential 0 , which is in turn the sum 
of external and induced potentials. The latter is simply 
taken as the Hartree potential produced by the perturbed 
electron density, while the former reduces to —r-E(t) for 
an incident electric field E(t). The induced dipole mo¬ 
ment is then calculated from the diagonal elements of 
the density matrix in the carbon-site representation as 
p(t) = —2e^i [pii{t) — where the factor of 2 ac¬ 

counts for spin degeneracy and = {xi^yi) runs over 
carbon sites. 

We use two different methods to solve Eq. and find 
p(t) (direct time-domain numerical integration and a 
perturbative approach m), which we find in excellent 
mutual agreement under low-intensity cw illumination. 
Direct time integration allows us to simulate the re¬ 
sponse to short light pulses for arbitrarily large inten¬ 
sity, while the perturbative method yields the nonlin¬ 
ear polarizabilities under multifrequency cw illumination 
(E(t) = Eq ee^^^^^ + Eo -\-c.c.^ with the same am¬ 

plitude Eq at both frequencies for simplicity), for which 
we can express the dipole moment as a power series in 
the electric field strength Eq according to 

n 

p{t)=Y^ E (So)” 

n si =—n 
S2 =—n 

+ C.C. (2) 

Here, n is the scattering order, whereas si and S 2 give 
the harmonic orders of the incident frequencies uJi and 
UJ 2 . We denote polarizabilities according to (^^’^^(co’out), 
where cjout is the frequency generated by a particular 
^th-order process. Considering terms up to third or¬ 
der (n < 3), Eq. [^then defines the linear polarizability 
{i = 1 or 2 ), the polarizabilities for SHG and 
THG, a^‘^\2uji) and (a^^^(3co’^), respectively, the wave¬ 
mixing polarizabilities corresponding to sum and dif¬ 
ference frequency generation ± ujj) {j = 1 or 

2 , i 7 ^ jf), and the four-wave mixing polarizabilities 
a^^\2uji ± ujj) and -h ujj — ujj). We obtain these 

polarizabilities by expanding the density matrix in Eq. 


as 

n 

p = E E pnSiS2^-i{siUJi+S20J2)t ^2) 

n si =—n 
S2 =—n 

which leads to self-consistent equations in (for the 

relevant combinations of n with harmonics Si and S 2 ). 
We solve Eq. following a procedure inspired in the 
random-phase approximation (RPA) formalism, as dis¬ 
cussed elsewhere [39] for SHG and THG. Eurther details 
on the extension of this formalism to cope with wave 
mixing are given in the Appendix. 

Acknowledgments 

This work has been supported in part by the Eu¬ 
ropean Gommission (Graphene Elagship GNEGT-IGT- 
604391 and EP7 -IGT-2013-613024-GRASP). 


Appendices 

The perturbative method used to simulate the optical 
response of nanographene to multifrequency continuous 
wave illumination is presented in detail. We also include 
additional results from numerical simulations performed 
in the time-domain, and provide further wave-mixing po¬ 
larizability spectra for graphene nanoislands with various 
sizes and dopings. 

Appendix A: Perturbative method applied to wave 
mixing 

Eollowing a procedure similar to that described in pre¬ 
vious work for the analysis of multiple-harmonic gener¬ 
ation [39], we express the single-particle density matrix 
equation of motion for a graphene nanoisland in the basis 
set of its tight-binding electronic states as 

=-'-^{^0 -£r)Pn' ( Al ) 

+ ^ E “ 4>i')ajiajiiipw - {pjj> - p %,), 

1 , 1 ' ^ 

where hsj is the energy of state \ipj) and = (/|0|/) 
are the matrix elements of the total electric potential in 
the basis set of the 2p orbitals |/) located at the carbon 
sites Rx The single-electron states and the carbon site 
orbitals are related through 

\(pj) = 10 ; (^ 2 ) 

I 

where the real-valued coefficients aji give the ampli¬ 
tude of orbitals |/) in states \(pj). These states are 
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FIG. 5: Linear absorption spectra for armchair nanotriangles. Spectra for triangular, armchair-edged graphene 
nanoislands with side lengths (a) L — 4.1 nm and (b) L — 8.4 nm are shown as a function of the number of additional charge 
carriers Q/e. The absorption spectra are normalized to the area of the triangle. We use a common color scale in (a) and (b). 


orthonormal QZidjidj'i = ^jj>) and form a complete 
set C^jdjidji' = thus facilitating transformations 
of the density matrix elements between site and state 
representations according to pjj' = '^ii> djidj'i'pw and 
Piv = ^jj' Pjj'' what follows, we use indices I 

to label carbon sites and j for singl e-electron states. 

In what follows, we solve Eq. ( |A1| ) for graphene nanois¬ 
lands exposed to multifrequency continuous-wave (cw) 
illumination described by the electric field 

E(t) = Eo + c.c.) e, (A3) 

where Eq is the field amplitude and e the polarization 
unit vector. Assuming that Eq is weak, we expand the 
density matrix as 


p — pnSiS2 ^-l{siOJi+S2UJ2)t 


(A4) 


j.nsiS2 

•jj' 


e 

h 


E E E 

n.' = ^ J —1 V 



SiUJi + S2OJ2 


where n = 1,2,3,... indicates the perturba tion order 
(i.e., terms proportional to (Eq)'^, see Eq. ( |A3| )) and 
Si (S 2 ) i s the harmonic index of oui (co’ 2 ). At 0^ order, 
Eq. (Al) is trivially satisfied with 
Eollowing Ref. [39], we insert Eq. (|A4|) into Eq. (Al) 


and collect terms with the same perturbation order and 
Q-i(siuji-\-s 2 UJ 2 )t dependence, from which we obtain the 
density matrix at order n > 1 as 


^nsiS2 

^jj' 


where 


® . „nsiS2 

h ^ SiWi + S 2 UJ 2 + ^1‘iT - (Sj - £j') 

(A5) 


in's'.s'o\ 

+ i/2T - [Sj — Ej') 


n—nbsi — si ,52 —S 2 

Pw ’ 


(A6) 


and 


f>T 


= - (R( ■ e) Eo5n,i (A7) 


i=l 


is the total potential. In Eq. (A7), the first term rep¬ 


resents the contribution from the external field (nonzero 
only for order n = 1), while the second term describes 


the Hartree potential produced by the perturbed electron 
density (here vw indicates the spatial dependence of the 
Coulomb interaction between electrons in orbitals |/) and 
|/')). The latter quantity ensures a linear dependence on 
pnsiS 2 fQj. term on the right-hand side of Eq. 


(A5), whereas lower perturbation orders are contained 


within • For each order n with harmonics si and 

52 satisfying n > |5i| + |52|, we are thus dealing with a 
self-consistent system in ^ which is treated using 
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FIG. 6: Plasmon-enhanced wave mixing of coincident pulses. We simulate the optical response of the nanoisland 
studied iu Fig. 1 of the maiu text {N = 1260 carbou atoms), for excitatiou by light pulses with the same duratiou (200fs 
FWHM) aud peak iuteusity (lO^^W/m^). Iu (a) we show the iuduced dipole momeut wheu the islaud is uudoped aud the 
iucideut pulses are tuued to the two most promiueut resouauces. Iu (b) we show the respouse for the same dopiug level 
cousidered iu Fig. 1 of the maiu text (three additioual electrous), but with cji aud CJ 2 tuued to differeut resouauces. 


the approach described in Ref. [39] . 

Specifically, we proceed by first using the identity 
Pw = J2jj' = Y^j O’jiO'ji'fj in the sum of Eq. 

( |A5| |, and then moving from state to site representation 
to obtain the diagonal density-matrix elements as 

E (A8) 

i' 

33 ' 


where 


Xw 


^ IP '' ■^^E + i/2r-(£,-e,,) 


(A9) 


is the noninteracting RPA susceptibility at frequency u. 

In summary, each new iteration order n is computed 
from the results of previous orders following a procedure 
similar to Ref. [39], here generalized to deal with multi¬ 
frequency cw illumination: 


1. We first calculate using Eq. ([A6). 


2. We then combine Eqs. (A7) and (A8) to find a self- 
consistent equation for , which reduces in ma¬ 

trix form to 




using site labels I as matrix indices and having de¬ 
fined 


— ~ (r^Z • g) Eo^n,l 


i=l 
nsiS2 


- 2e vwajiajniljp 
I'jj' 


3. We use the calculated values of and 

to obtain using Eq. (A8), and from here the 


induced charge at site I at order n associated with 
the harmonics si and S 2 as Com¬ 

putation and storage demand are reduced by using 

the property ' 


4. Einally, the polarizability for wave mixing among 
the harmonics si and S 2 is calculated from 


{siUJi + S2U32) = 


2e 


{EoY 


(AlO) 


upon iteration of this procedure up to order n. 


Extension to more frequencies is straightforward. 
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FIG. 7: Difference frequency generation. The nonlinear polarizability —UJ 2 ), corresponding to difference frequency 

generation, is shown for a triangular, armchair-edged graphene nanoisland containing N = 330 carbon atoms and doped with 
(a) three and (b) six additional electrons. In each case the upper panels in (a) and (b) show the linear absorption cross-section 
of the nanoisland, while the lower panels show the polarizabilities as the incident held frequencies ui and UJ 2 are varied. 
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FIG. 8: Four-wave mixing. The nonlinear polarizability a^^\2uji + 002 )^ corresponding to four-wave mixing, is shown for 
the same nanoisland and doping levels considered in Fig.|^ 


Appendix B: Effect of doping on wave mixing 

The linear absorption spectrum of an armchair-edged, 
doped graphene nanoisland exhibits numerous plasmonic 
resonance features, some of which display strong electri¬ 
cal tunability. In particular, the lowest-energy prominent 
peak, which is absent when the nanoisland is undoped, 


undergoes dramatic energy shifts upon adding only a few 
additional charge carriers (see Fig. |^. As this mode is 
tunable, and also tends to produce large nonlinearities 
m, we excite it in the nanoisland considered in Fig. 
1 of the main text with one of two collinear ultrashort 
pulses to study wave mixing. In Fig. we investigate 
wave mixing of collinear pulses in the same nanoisland, 
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but here one of the pulses is tuned to a different plas- 
mon from that of the main text. Interestingly, a different 
qualitative behavior is observed when the island is ei¬ 
ther undoped (Fig. [^) or doped (Fig. |^). The former 
indicates that the second-order response vanishes for po¬ 
larizations both parallel (upper panel) and perpendicu¬ 
lar (lower panel) to an edge of the nanoisland when it 
is undoped, while a strong third-order response remains, 
although it is weaker than that of Fig. lb. When the 
nanoisland is doped, the second-order response is recov¬ 
ered for the perpendicular polarization, but overall the 
nonlinear response is weaker than that of Fig. lb, where 
the highly-tunable mode participates in wave mixing. 

Appendix C: Additional wave mixing polarizabilities 

In Figs. and we present the wave-mixing polar¬ 
izabilities corresponding to difference frequency gener¬ 
ation, a{uji — UJ 2 ) and four-wave mixing, a{2uji + co’ 2 ), 


respectively, for the N = 330 graphene nanoisland con¬ 
sidered in Figs. 2 and 3 of the main text doped with 
(a) three or (b) six electrons. As before, we find sim¬ 
ilar plasmonic enhancement of the input frequencies in 
both cases (see horizontal and vertical features in the 
contour plots), while in Fig. enhancement at the out¬ 
put frequency follows the curve Co ’2 = Co’i Tco’p and in Fig. 
I^the output enhancement follows co ’2 = — 2co’i + co’p, uOp 
being any of the plasmonic resonances in an island. In 
Fig. 13 we note a very large difference frequency genera¬ 
tion polarizability near huui ~ 1 eV, huj 2 ~ 2 eV (or vice 
versa), where in fact we have a triple-resonance condi¬ 
tion: simultaneously, both of the input frequencies and 
the output frequency are resonant with plasmons in the 
nanoisland. In Fig. the cases where co’i = co ’2 = 
correspond to plasmonic enhancement of third-harmonic 
generation (THG). For comparison, the third-harmonic 
susceptibility (3co’) has been measured in 10 nm silver 
nanoparticles as ^ 2 x 10“^^ esu [56], corresponding to a 
third-harmonic polarizability (3co’) of ^ 2 X 10 esu. 
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